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Abstract. We study classical binary fluid mixtures in which densities vary on very 
short time (ps) and length (nm) scales, such that hydrodynamics does not apply. In a 
pure fluid with a localized heat pulse the breakdown of hydrodynamics was overcome 
using Green's functions which connect the initial densities to those at later times. 
Numerically it appeared that for long times the results from the Green's functions 
would approach hydrodynamics. In this paper we extend the Green's functions 
theory to binary mixtures. For the case of isothermal isobaric mutual diffusion in 
isotopic binary mixtures and ideal binary mixtures, which is easier to handle than 
heat conduction yet still non-trivial, we show analytically that in the Green's function 
approach one recovers hydrodynamic behaviour at long time scales provided the system 
reaches local equilibrium at long times. This is a first step toward giving the Green's 
function theory a firmer basis because it can for this case be considered as an extension 
of hydrodynamics. 
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1. Introduction 



With the increase in miniaturization nanoscale technology is developing rapidly. 
However, at such small scales the physics of thermodynamic properties and transport 
processes will deviate from that at larger scales. Certainly one sees finite size effects and 
appreciable fluctuations on measurable quantities. It is likely that new physics exists at 
the nanoscale. 

Of course, on large time and length scales transport processes in fluids are well 
described by hydrodynamics which contain the transport coeflicients such as the 
diffusion constants, viscosity and thermal conductivity [1]. Kinetic theory can provide 
low density approximations to these coefficients [2] while linear response theory yields 
the Green-Kubo relations which express the transport coefficients in terms of integrals of 
time-correlation functions of the appropriate current ffuctuations in equilibrium [3,4]. 
Thus once one has determined the transport coefficients by studying the equilibrium 
time-correlation functions, one can attack a wide range of non-equilibrium transport 
problems described by hydrodynamics not too far from equilibrium. In contrast, the 
description of the non-equilibrium behaviour at shorter times and smaller length scales 
remains mostly unsolved despite its importance, in principle, to the understanding of 
nanostructures. Even low density kinetic theories cannot describe time scales below 
the average duration of a collision. Attempts to generalize hydrodynamics stranded 
at first because a straightforward further expansion in the gradients, i.e. beyond the 
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linear laws, leads to higher order Burnett coefficients which are infinite [5]. Recently 
a theory was proposed by Prof. J.M. Kincaid using Green's functions to describe 
processes at small length and time scales which avoids these divergences [6-9]. Like 
hydrodynamics the Green's function approach can also treat the hydrodynamic fields 
(density, momentum and energy) but unlike in hydrodynamics currents need not be 
hnear in the gradients. Another contrast with hydrodynamics is that the Green's 
function approach is not restricted to long times. Thus the Green's function approach 
should be able to attack a wider range of problems than standard hydrodynamics. 
Indeed the Green's function approach has been numerically apphed with success to self- 
diffusion [6] and heat transport [7-9] on nanometre and picosecond scales. In these 
applications, the Green's functions are resummed into an infinite series that needs to be 
truncated. It was shown using computer simulations that if just a few terms from the 
infinite sum are taken into account, accurate results were obtained for all time scales for 
those systems [6-9]. Also, for long times the numerical results were seen to approach 
the numerical hydrodynamic results [9]. 

In this paper we formulate the Green's function approach for mixtures, focusing on 
mutual diffusion in isotopic and ideal binary mixtures. Apart from adapting the Green's 
function approach for mass transport in mixtures, this paper is mainly dedicated to the 
question whether the truncated expressions for the Green's functions describe the right 
very short time behaviour and whether they do indeed reduce to hydrodynamics at 
long time scales and for small gradients; if the latter is true then the Green's function 
approach can rightfully be called an extension of hydrodynamics. For self-diffusion in a 
single component fiuid the truncation was already shown analytically to give the right 
hydrodynamic- like behaviour for long times [6] . Mutual diffusion in an isotopic or ideal 
mixture is the simplest case to consider next (simplest both in thermodynamic and 
in hydrodynamic aspects) Ideally one would like to show analytically that the Green's 
functions reproduce not just mutual diffusion but the full hydrodynamic equations in 
all cases where gradients are small and times arc long, but this would be an enormous 
undertaking. Therefore, as a first step, we restrict ourselves in this paper to the already 
nontrivial case of diffusion in isotopic and ideal mixtures. 

2. Classical isotopic and ideal binary mixtures 
2.1. Definitions 

We consider a classical mixture of two components where the particles interact through 
short-range pair potentials. The total system consists of N particles of which A^i 
particles are of component 1 and N2 — N — Ni particles are of component 2. The system 
is enclosed in a volume V. We denote the (three-dimensional) position of particle i of 
component A as r\i and its velocity as V\i, where the index A equals 1 or 2, and i is an 
index which runs from 1 to Nx- The mass of particles of component A will be denoted 
as mx- We restrict ourselves to two related classes of mixtures: isotopic mixtures, for 
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which the pair potentials are the same independent of the components to which the 
interacting particles belong, and ideal mixtures, whose pair potentials can be neglected 
as far as thermodynamic properties are concerned. 

The system is prepared in a non-equilibrium state described by some given, non- 
equilibrium initial ensemble in phase space 'Pne(r), where T = {r\i, V\i}. In the initial 
non-equilibrium ensemble — and consequently also at later times — the particles of each 
component are indistinguishable. After the initial preparation in the non-equilibrium 
state, the system stays isolated, i.e. the system does not interact with the outside. 

The macroscopic (i.e. averaged) local density of component A is defined as [10] 

n^{r,t) = / drVU^)J2^ir - rxiit)) = N^{S{r - ni{t))U (1) 

i=i 

where the subscript "ne" on the average indicates its non-equilibrium nature. In this 
equation we used the indistinguishability of the particles of the same component to 
express the densities nx{r,t) in terms of the properties of a representative particle of 
each component i.e. in terms of rxi{t). 

Due to the non-equilibrium nature of V^e, the local densities ^^(t', t) will in general 
not be uniform at t = 0. In this paper we are mainly interested in the time evolution 
of these densities. This is in principle determined by the microscopic time evolution 
governed by the Hamiltonian 

2 Nx 

^(n = EE>Ait^A.r+^(r) (2) 

A=l i=l 

where 

C/(r)= 5^'i^,,(|r,,-r.,|). (3) 

X,u,i,j 

Here (fx^, are the pair potentials and the summation ^ is restricted to exclude self- 
interaction of the particles, i..e. to exclude the cases where both X — u and i — j. For 
isotopic mixtures, (pu — (fin — (fi22- The equations of motion are thus given by 

rxi = vxi (4) 

^Ai = ^ — ■ 5) 

mx orxi 

Before proceeding we introduce a few more quantities for later convenience. The 
local mass densities are given by px ~ rnxux- Here and below when the spatial and 
temporal arguments of a local quantity are not denoted, it is implied that they are taken 
at position r and time t. The local total number and mass densities are, respectively, 
n = ni+n2 and p = Pi+P2- We furthermore define the local mole fractions as Xx = nx/n 
and the mass fraction of component 1 as c = pi/p. The mass fraction of component 2 
is simply 1 — c. The local mass currents of the two components in the laboratory frame 
are given by (A = 1, 2) 

fx = Nxmxivxiimr - rAi(t)))ne (6) 
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Using these definitions the local mass velocity of the fluid can be written as 

^=(ji*+J2*)/P- (7) 

Finally, we define the diffusive currents of the two components in the local frame (i.e. 
co-moving with the local velocity u) as 

JA = J A - Pxu = Nxmx{[vxi{t) -u]6{r - rxi{t)))ne (8) 

so that by definition we have 

ji+j2 = 0. (9) 

2.2. Thermodynamics, hydrodynamics and Green- Kubo formula 

For classical isotopic mixtures and ideal mixtures, many thermodynamic equilibrium 
properties take on a particularly simple form — which is in fact the main reason we 
consider these systems here. In particular, for the current case of mass transport 
investigated in this paper, we will need {dfi/dc)Tp and {dp/dc)Tp where T is the 
temperature, p is the pressure and 

/i = (10) 

mi 1712 

and fix is the chemical potential of component A. These thermodynamic properties can 
be found from the form of chemical potential for isotopic systems [11, §93]: 

Hx = kBT\ogxx + fiQi= fio2- ^ksT log — , (11) 

7712 

where /Xqa is the chemical potential of the pure isotope A at the same temperature T 
and pressure p as those for which fix is evaluated. Equation (|TT|l implies that 

'dfix\ ksT fdfii\ fdfi2 



OxxJtp XX ' \dp Jtx Jtx' ^^^^ 

We can use the Gibbs-Duhem relation [12, §8.7] 

Nidfii + N2dfi2 = Vdp - SdT (13) 

to derive that Ni{dfii/dp)Tx + N2{dfi2/dp)Tx = V- Combined with {dfii/dp)Tx = 
{dii2/ dp)Tx (cf. equation (fT^ ) this yields 

We can write c as a function of Xi as c = miXi/[miXi + m2(l — Xi)], which gives 
dxi/dc = p^/(mim2n^). Using also equation (|Tn|l . we get for equations (|T^ and (|T^ in 
terms of c and /i 

'dfi\ _dxi/dji\ _ pksT 
dcjj,^ do V^^i/Tp c{l — c)mim2n 
dfi\ f 9fi\ m2 — mi 



op/Tc \9pJtx niim2n 
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From the fundamental relation for the Gibbs free energy dG = —SdT+ Vdp + pidNi + 
P2dN2 [12, §8.7], one can deduce that {dp/dc)Tp = —p'^{dfi/dp)Tc, whence 

'dp\ ^ (mi - m2)p^ 
dcjj,p mim2n 

The quantities in equations (P3|) and (fTTjl will be needed in the treatment of the 
long time behaviour of the Green's functions for mass transport in section |3J These are 
precisely of the form valid for ideal ("perfect") mixtures [11, §92], which is the reason 
why our treatment is valid both for ideal mixtures and for isotopic binary mixtures. 

Because in section H] the results of the Green's function approach for long times 
will be compared to those of hydrodynamics, we will now review the appropriate 
hydrodynamic equations and Green-Kubo formula for mass transport in ideal mixtures 
and in an isotopic binary mixture. Since we restrict ourselves in this paper to mass 
transport, the relevant conservation laws are only those for the local densities n^, which 
follow from equations ((H), (jH)) and (jHl): 

^ = - V • 3l/mx = -V ■ {nxu + jx/mx) (18) 
at 

In terms of p = mirii + m2n2 and c = mini/p, these equations become [1, ch. 2] 

|=-PV.„ (19) 

pg=-V.i, (20) 

where we have used the material derivative D/Dt = d/dt + u ■ 'V. To make these 
equations into a closed set, ji needs to be expressed as a function of p and c (as well as, 
in general, of u and T, which will be assumed to be constant below). Near equilibrium 
one can use linear laws as closure relations. Following a widely used choice of currents 
and forces in the theory of irreversible thermodynamics [1, ch. 3], one has in general 

ji= - ^[Vp]t-^V\^T (21) 

Here the gradient [Vp]t is taken at constant temperature while the pressure p and the 
concentration c can vary. Ln and Liq are called phenomenological coefficients. 

We will only consider here the case in which there is only a gradient in the 
concentration c, while gradients in temperature T, pressure p and fluid velocity u are 
zero. In that case, one can write equation (j2T| as 



L 



Here the mutual diffusion constant D = {dp/dc)TpLii/ (pT), or, with equation (|T5|l for 
binary isotopic mixtures and ideal mixtures, 

D = . (23) 

c(l — c)mim2n 



dp\ 
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dt' (MO) ■ Mt')) (24) 



The Green-Kubo formula for the phenomenological coefficient Ln for mutual 
diffusion in a fluid is [3,4]: 

where the average is over a grand canonical ensemble in the thermodynamic limit and J\ 
is the microscopic expression for the total current in the centre of mass frame given by 

JA(t)=mA5^[«A.(t)-*] (25) 

i=l 

in which v is the centre of mass velocity, i.e. 

2 Nx 
A=l i=l 

The Green-Kubo formulae for Liq and other transport coefficients are more complicated 
and not needed here; for these we refer to the literature [3,4]. 

Note that v is not an ensemble average but an average over particles taken per 
ensemble member, and that it is a constant of the motion for that ensemble member. 
With the inclusion of v we have on a microscopic level the analogue of equation i.e. 
Ji + J2 = 0. This relation can be used to find an alternative representation of Lu: 



^11= / dt' (MO) ■ J2{t')) (27) 

3k bV Jq 

Using equations and ^I7\i in equation for D, and that particles of the same 
component are indistinguishable, we flnd that 

1 f°° 

^ = ~ 3c(l - c)nV Jo - ^] ■ [^21 - (28) 

In the thermodynamic limit (A'^i, N2, V ^ 00 keeping Ni/V and N2/V flxed and flnite) 
we may factor the averages to leading order, so that we get 

D = ^/i2, (29) 

where [13] 

1 r°° 

fi2 = 3 dt' (iV[«n(0) - v] ■ [V2i{t') - v]). (30) 
Similarly, for future use, one can also deflne, following reference [13] 

fxap = 3 dt' (iV[t;Ai(0) - v] ■ [vx2{t') - v]) (31) 

Dx =1 r^t' {[vxiiO) - v] ■ [v,,{t') - v]). (32) 



3 _„ 

Note that the subscript a/3 does not denote indices but instead signifles that correlations 
between different particles are involved and that we added a factor N in equations (|3U|) 
and compared to the deflnition of /12 and fxa/3 in reference [13] to get a well-deflned 
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thermodynamic limit. Due to the subtraction terms —v and the conservation of total 
momentum, one has the relations [13] 

fi2= - -^{Di + xih^p) (33) 

m2X2 

fl2= --^{D2 + X2f2ap) (34) 

miXi 

Using these, one can rewrite the right hand side of equation (pUj) as [13-15] 

D = X2D1 + X1D2 + XiX2(/la/3 + /2a/3 - 2/12). (35) 

The first part of this relation, i.e. the part involving only the self or intrinsic diffusion 
constants Dx, is often used as if it were the full expression for D, which goes under 
the name of the Darken [16] or Hartley-Crank equation [17]. Douglass and co-workers 
were apparently the first to have shown that corrections to the Darken-Hartley-Crank 
equation involve the cross correlations given by the fxa/3 and /12 [13,14]. 

3. Green's functions for mass transport in mixtures 
3.1. General formulation 

The hydrodynamic treatment only works near (local) equilibrium. For other situations 
i.e. on small length and time scales, a new description is necessary. For self-diffusion 
the short time behaviour, as well as the long time behaviour, could be successfully 
described by using Green's functions [6]. A similar description for heat transport gives 
good numerical agreement with simulations as well [7-9]. 

Analogous to the Green's function approach for self-diffusion we want to write the 
r and t dependent densities ux here as integrals over Green's functions that give the 
probability for a representative particle of the appropriate component to be at position 
r at time t given that it was at position r' at time 0. Starting from equation (^J) we can 
write (A = 1,2) 



nx{r,t)=Nx{6ir-rximnc= / dr' Gx{r,r' ,t) nx{r' ,0), (36) 

Jv 

where the Green's function Gx is defined as 

Gx{r, r', t) = {6{r - r' - ^rxi{t)))xr' ■ (37) 

Here ArAi(t) = rxiif) — r'Ai(O) is the displacement of particle 1 of component A while 
Oaw is the average over initial conditions for which t'ai(O) = r'. For an arbitrary phase 
function A(r), this conditional average is defined as 

(„r,),., . <-^'^'^"-'-7;'°»>- . (38, 

(d(r' - rAl(0)))ne 

We now make two remarks about this definition of the Green's functions and the 
conditional average. First, for initial conditions where nxir' , 0) = A^a(^(^'' — ■'"Ai(0)))ne = 
for some r', the right-hand side of equation ()38|) seems ill defined. The way we choose 
to define it for such r' will not matter since (5(r — r' — Arxi{t)))xr' is multiplied by 
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nx{r',0) = in equation anyway.l Second, equation differs in normalization 
from the definition of ( ),.' in the heat pulse case of references [7-9], in which the 
conditional averages were only properly normalized when the system is initially of 
uniform density (which was the case in references [7-9]). In contrast, the definition 
of the conditional average in equation ()38j] is always normalized i.e. (l)Ar' = 1- This 
definition is also more analogous to the self-diffusion case [6] because of the appearance 
of nx{r', 0) inside the integrals in equation (^^1 instead of a single factor of the average 
density N\/V in front of the integral. 



3.2. Expansion around Gaussian behaviour 

Next the Green's functions Gx are Fourier transformed with respect to r — r', i.e. one 
considers 

Fx{k,r\t)= [ dre"'<^-^'^Gxir,r',t) (39) 
Jv 

Using equation and expanding in k, this gives 

Fxik,r',t) = {exp[ik ■ Arxiimxr' = V ^yUA„(fc, r', t) (40) 

n=0 

where k = \k\, k = k/k and the moments are defined as 

fiXnik,r',t) = (\k-Arxiit)Y) . (41) 

\ L J / Ar' 

Since divergences would occur if one would directly take the inverse Fourier transform 
term by term in equation (^UJ, an alternative expansion will be used. For given direction 
fc, Fa is the generating function of the cumulants of the displacement of the particle in 
that direction, [18-20] so that one can write 

Fx{k, r', t) = exp V ^^Kxn{k r\ t) (42) 

n=l 

The cumulants Kxn and ^ixn are related through [20] 

{p^>0} 1=1 1=1 ^ 

where we dropped the arguments fc, r' and t of ^xi and Kxn, as we will do below as well. 
For instance, for the first few ^An, equation ()43p becomes 

«A1 = yUAl (44) 

K.\2 = /iA2 - fJ'll (45) 

f^X3 = /^AS - 3/iAl/iA2 + 2/1^1 (46) 

KA4 = /iA4 - 4yUAlyUA4 " 3/i^2 + 12/iA2/UAl " ^fJ'XV (47) 

if However, to have a proper definition everywhere, one could follow the convention that {A)\r' = {A)na 
for values of r' where nx{r',0) — 0. An alternative would be to introduce a "tracer particle" of 
component A at r' but this works only in the thermodynamic limit. 
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The Fourier transform in equation (j42j) is first rearranged by factoring out a 
Gaussian expression [6-9]. This can be done in several ways, but here we write 



F,{k,r',t) = Mk,r',t) 



n=3 



(4^ 



and we demand that the first correction term to the Gaussian f\{k, r\ t) in equation 
is only 0{k^), so that 

fx{k, r', t) = exp [iKxik - ^Kx2k^] (49) 

with Kxi the mean of the displacement and k\2 its variance, explicitly given by 

Kxi = {k- Arxi{t))xr' (50) 
Kx2 = {[k- Arxiit) - Kxiik r\ t)]^)xr'. (51) 

The coefficients bxn in equation ()48|1 can be determined by gathering terms of like order 
in k in the fc-expansion of Fx{k, r', t) / fx{k, r', t). This gives 



E n 

{Pi>0} ^=3 



' 1 









(52) 



where the arguments k, r' and t of bxn and Kxn have been suppressed. The first few 
coefficients bxn are given by 

bx3 = «:a3/3! (53) 
&A4 = «:a4/4! (54) 

Up to a constant prefactor, the coefficients kxi (mean), kx2 (variance), bx3 (skewness), 
bx4 (kurtosis) and bx5 are pure cumulants, however starting from n = 6 the coefficients 
bxn can only be written as combinations of cumulants, e.g. bxe = {kxq + 10k|3)/6!. 

Term by term the expression in equation (j48j) can be Fourier inverted. However, 
one has to keep the k dependence of the Kxn in mind. For this k dependence, we can 
write, using equations (HH), (^Hj) and (j^ . 



Z^An 
^An 
bxn 



KAn 
BAn 



(55) 
(56) 
(57) 



where MAn, KAn and BAn are r' and t dependent tensors of rank n. Here a tensor 
notation is used. In this notation, ":" denotes a scalar product (summation over all n 
indices), while powers of tensors are A" = A A - ■ ■ ® A and yield symmetric tensors 
of rank n. The tensors Mxn, ^xn and BAn are chosen to be symmetric, and follow from 
equations (gH), (0, (US) and equations (IHHll- (IK7|l :^ 

MAn= ([ArAl(t)]")Ar.' (58) 

§ To avoid having to use explicit indices, we use here a symmetrized tensor product, such that 
AB = [A (g) B]s, where []s denotes the symmetric part of a tensor e.g. [A]sy = {A\j + Aji)/2! for 



a rank two tensor, [AJsyfe = (AajTc + A 



■jik 



l^kij + Ajfei)/3! for rank three, etc. 
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K 



An 



{pe>0} ^=1 ^=1 



Ma,/£!; 



B 



An 



K 



A£ 



E n 

{pe>0} £=3 

T.T=-i^Pi='^ 

Using these definitions, the Fourier inverse of f\ of equation ()49p is 

exp[-i(r - r' - Kai) ■ K;^2' ■ {r - r' - Kai)] 



(59) 



(60) 



9x{r,r',t) 



Vdet(27rKA2) 



(61) 



where 



Kai = (ArAi(t))Ar' (62) 

KA2 = ([ArAl(t)-KAl]')Ar'. (63) 

Using equations (jiH|) and (jHTj) and that the Fourier inverse of Fa is Gx and that of fx is 
gx, we obtain 



Gx{r,r',t) 



1 + ^ Ba„ 



n=3 

oo 



1 + ^BAn: H„(r-r'-KAi, ^K-^) 



n=3 



(64) 



gxir,r',t). (65) 



Here the H„ are polynomial tensors in s = r — r' — Kai given (with A = |K^2^) by 



H„(5,A) = ^ 



-l)™n! 



m=0 



m\{n — 2m)! 



A'"(2A ■ s 



in— 2m 



(66) 



where [^J is the largest integer number less than or equal to ^. Explicitly, the first few 
polynomials H„>3(s, A) are given by 

H3(s,A) =4[2(A-s)=^-3A(A-s)] (67) 
H4(s, A) = 4[4(A ■ - 12A(A ■ + 3A^]. (68) 

The coefficients in the expressions of the right hand sides for these H„ are precisely those 
of the Hermite polynomials. They can in fact be related to the Hermite polynomials, 
in that the nth Hermite polynomial can be expressed as i7„(s^) = H„(s, 1) : x"^ (using 
Gradshteyn and Ryzhik 8.950.2 [21] with s = [s^, Sy, s^), 1 the identity matrix, A : 
the (all) X component of A, and noting that 2'"(2m — l)!!/(2m)! = l/m\). 



3.3. Reduction to one- dimensional Green's functions 

The somewhat cumbersome tensorial notation of the last section can be avoided in 
cases where the initial distribution Vne{^) is only nonuniform in the x direction and 
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translation invariant in the y and z directions, such as in previous work on the Green's 
functions [6-9]. In such a situation one can use the reduced Green's function 



X, X , t) 



jdy'jdz' Gx{r,r',t) 



{6{x -x' - Axxi{t)))xx' 



(69) 



where Axxiit) = xxi(t) —xxi{0) is the displacement of particle 1 of component A in the x 
direction. Note that the y and z dependence has dropped out because of the translation 
invariance of Pne(r) in those directions and the translation invariance of the dynamics 
i.e. of the Hamiltonian. In terms of the reduced Green's function equation becomes 

nx{x,t)= j dx' Gl{x,x',t)nx{x',0). (70) 

To see how the reduction works out for expansion (jHS)), it is easier to start from the 
original expansion in equation ()64j). Noting that because of translation invariance the 
Kxi{r',t) and consequently the Bxni'f',t) do not depend on y' or z', we can integrate 
equation (jMj) over y' and z' to find 



where 



Bxn '■ x"- and 



oo ^ 



d_ 

dx 



gl{x,x\t) 



(71) 



gl{x, x\ t)= dy' / dz' gx{r, r', t) 



exp 



(X 



Kxn '■ and we used that the xy and xz components of the matrix K;v2 
are zero by translation symmetry in the y and z direction (which also pertains to the 
inverse of this matrix). Working out the derivatives to x in equation ()7H) gives 



where k 



Xn 



X — X 



K 



Al 



glix,x',t), 



(73) 



„,=3 (2«^A2)"/' 

where the Hn are true Hermite polynomials. 

We note that to obtain the and k,^^ it is possible to avoid the tensors through 
which they have been defined, since by comparing the definitions = Kxn '■ and 

= Bxn '■ with the expressions in equation (jKHjl and equation (jKTjl . one sees 
immediately that one only needs to set k equal to x in the expressions for bxn and Kxn 
to obtain b^^ and 

We note that the resummation of Gx leading to equations and (fTSj) is somewhat 
different from that used for the non-equilibrium heat transport case of references [7-9] 
where separate Gaussian factors were introduced for the even terms and the odd terms 
in k (i.e. even and odd n), respectively, rather than a single overall Gaussian factor as we 



used in equation ()48|1 . A result was then found in terms of Sonine polynomials S", 



X 



[2]. 



The present formalism is nonetheless very similar since the Hermite polynomials are 
special cases of the Sonine (or Laguerre) poljTiomials (i.e. H2n{x) = (— l)"2^"r7,!S'i"''y2(^^) 
and H2n+i{x) = (-l)"22"+in!a;5{/^(a;2), see references [2, p. 127] and [21, 8.972]). 
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Leaving out the odd terms in equation (|73p furthermore yields the exact same form 
as that in the self-diffusion result in equation (6) of reference [6] (with pi = ^1^x2^ 
P2 = 6^4 and p3 = 5^6)- 

3.4- Practical considerations 

Substituting the expression for Gx in equation into expression (jHUj) (or equation (f75|) 
into equation (HO))) for the local densities nx in terms of the Green's functions, leads to 
an expression for those densities involving an infinite sum over n. To make the Green's 
function approach practical the dynamical averages occurring in Mxn in equation 
(and in Ka„ and Bxn through equations (jSUj) and (jUUj) ) need to be calculated, for instance 
by molecular dynamics. Since one can calculate only finitely many quantities in practise 
the infinite sum in equation ()65|) has to be truncated at some point. 

Having chosen the Gaussian factor gx in equation ()48|) such that corrections are at 
most 0{k^), the Gaussian approximation 

Gxir,r ,t) ^ gx{r,r ,t) = .(74) 

^det(27rKA2) 

(cf. equation (jblj) ) is expected to describe the correct long time limiting behaviour, 
similar to the self-diffusion case [6] . To confirm this expectation we will look in detail at 
the behaviour in the long time limit in section |31 For intermediate (non-infinitesimally 
small) times, it will be necessary to include some of the correction terms proportional 
to Bxn in equation (j65|l e.g. (using equations (j58|l - (j60|) . (j67j) and 



Gxir, r', t) ^ gx{r, r', t) [l + K,3 : {ql - ^^iM 

+ Kx,:{q{-QK-xlql + ?>[K-xlf)+ ...], (75) 

where = ' — ~ ^ai). Checking whether these first few terms for the Green's 
functions describe the intermediate time behaviour accurately will be done in future 
numerical work. 



3.5. Connection with the Van Hove self- correlation function 

For equilibrium initial conditions (i.e. replacing ()xr' by ()), the Green's function 
Gx{r,r',t) depends only on the difference \r — r'\ and reduces to the Van Hove self- 
correlation function [22,23]. Indeed, the self-part Gs{r — r', t) of the Van Hove function 
also has the physical interpretation of the probability for a particle to be at position r 
given that it was at position r' a time t earlier. Thus the Green's function Gx{r,r',t) 
can be seen as a nonequilibrium extension of Gs{r — r',t). It is at present not clear 
whether this means that Green's functions could be observed in incoherent neutron 
scattering experiments like the Van Hove self-correlation function can [22]. 

Similar expansions as in equation ()48j) have been used for the Van Hove self- 
correlation function Gg in a single component fluid in equilibrium [20,23-26], where odd 
terms in k are absent. The Gaussian approximation to the Van Hove self-correlation 
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function turns out a good description in the very short and the long time limits. 
Non-Gaussian corrections are important for intermediate times, where including just 
the first correction term (involving the fourth moment) captures most of the physical 
behaviour [23-25]. For more information, we refer to chapter 4 of reference [23] which 
contains a nice overview of different approaches to Gg- 

4. Limiting behaviour of the Green's functions at long times 

As a first approximation in investigating the large time and length scale behaviour of the 
Green's functions, we realize that at long times gradients get small so that 0{k^) terms 
become unimportant and the Gaussian approximation, found by neglecting the terms 
proportional to Bxn in equation {0{k^) terms in equation (jIH|) ) should capture all 
of the relevant dynamics for long times, since the hydrodynamic equations also contain 
at most second orders in gradients. This will now be checked. 

Taking into account only the Gaussian term in the resummation in equation 
we get the Gaussian approximation in (fflj) . with Kai and Ka2 given in equations (jH^ 
and ()63p . We emphasize that the present approach differs from the self-diffusion case [6] 
where no odd terms in k existed. But in the mutual diffusion case odd terms such 
as KAi(r', t) cannot be neglected. These terms are of the form {Arxi(t)) xr' so they 
represents net particle displacement or drift. In a diffusing system particles do move in 
a preferred direction, namely opposite to the gradient, so these terms will not be zero 
and will in fact be sensitive to the gradients. 

Using the Gaussian approximation Gx ~ gx in equation (jHj), a diffusion-like 
equation can be derived for the Green's function Gx- The easiest way to derive this is 
by taking the time derivative on both sides of equation P9|l . thus obtaining an equation 
for the Fourier transform Fx ^ fx, and then taking the inverse Fourier transform to get 
the corresponding equation for Gx ^ Qx- This leads to 



where Kai = ^^Kai 
this implies that 

dnx 




at 



\ [ dr'kx,ir',t)Gxir,r',t)nxir',0) 



Jv 




(77) 




(78) 



Jo 
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(using the symmetrized tensor product explained in the footnote on page II U|). Equation 
(f7S|) looks surprisingly like a Green-Kubo formula for the self-diffusion coefficients except 
for the conditional average {)\r' and the integration to t. It seems reasonable to assume 
that for long times t the condition that a particle was at some position r' at time 
should have a negligible effect on how fast it is diffusing at time t. In other words it 
seems reasonable that the limit of t — > cxd of the above expression is independent of r': 

Mm k^2{r',t) = hm 2 f dt' {[v^,{t) - («Ai(t))][«Ai(t') - {vxi{t'))]) (79) 

For if: — ^ oo the (closed) system will have reached equilibrium, so we can use the time 
translation invariance and isotropy of the equilibrium distribution and that {v\i{t)) = 
for all t, to write 

2 r°° 

limKA2(r',t) = - / dt' {vxiiO) ■ vxiit'))l = 20^1, (80) 

where we used the definition of D\ in equation ((221) • (Note that formally one needs to 
take the thermodynamic limit before taking t oo.) 

Assuming now that Kx2{r',t) attains its limiting value 2Dx fast enough, we can 
replace [<.\2{r',t) by that value in equation (f77|) for long times. It is then independent 
of r' and can be taken outside of the integral in equation ((77j). (Another way to see 
that the r' dependence of K\2 should drop out is by noting that this r' dependence 
will just give additional gradients in r if one expands around r, but the Ka2 term in 
equation (f77|) is already of the highest order in gradients required for hydrodynamics, 
i.e. second order in the gradients.) Equation (|77|l then becomes 

- V ■ J^dr'k»ir',t)Gxir,r',t)n,ir',0) + D,VW (81) 

Because the first expression on the right-hand side of equation ()8H1 cannot be 
further calculated without additional assumptions, to see whether it describes long time 
behaviour in accordance with hydrodynamics, one next makes the essential assumption 
of local equilibrium (which is also used to derive hydrodynamics in other ways). 

To keep the derivation of hydrodynamics from the Green's functions simple we will 
also assume that the system is in mechanical equilibrium i.e. Vp = 0, and only has 
gradients in the concentration i.e. VT = and Vn = 0. 

Within the context of local equilibrium, assumed to hold at time zero for now, 
we focus on the quantity K;vi in equation (|HT|) for component A = 1. According to 
equation ((U^ . 

''^ = {Sir' -mm).. = M^) ' ^''^ 

where we have used that particles of the same component are indistinguishable. We will 
now evaluate Kii(r', t) with the local equilibrium distribution Pl derived in appendix A 
(following Ernst's treatment of a single component fluid [27]). Using equation ()A.13|) 
we can write Kii(r',t) as 
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Kii(r',t) 

' N,v^,{t)6{r' - rn(0))|l + p [ dr" l^^l^ _ J'^^^'P 1 (^r" -r + ut)-Vc] 
\ I y cm2n (1 — cjrriin J 



ni(r',0) 

(83) 

where (yl)o is defined as / A(T)Vo{T)dT where Vo is the grand canonical distribution but 
with all velocities shifted over u, and nx{r") is the microscopic expression for the number 
density of component A at position r" (see appendix A for details). In appendix B it is 
shown that this gives to leading order in the gradients: 



(iVi(iVi - l)t;n(t)ri2(0)) (iViiV2^n(t)r2i(0)) 



■ Vc. (84) 



cm2nV (1 — c)minV 

Using now that (i'ii(t)ri2(0)) = — (i'ii(0)ri2(t)) = dt' {vii{Q)vi2{t')) and furthermore 
taking the thermodynamic limit and then the limit t — >■ oo, we get (cf. equations (|3U|) 

____ 2 

and dnH)) Kii(r',t) =u + ^^^^„2 (/12 - /iq/3)Vc. By interchanging 1 and 2, which also 
imphes that Vc be replaced by V(l — c) = — Vc, we find the corresponding result for 

2 

component 2, i.e. K2i(r',t) =u - ^^^^„2 (/i2 - /2a/3)Vc, so that 

kAi(r',t) = n - (-1)" , (/i2 - /a„/3)Vc (85) 

Substituting this result for Kxi from equation (|H3j) into equation (jHT|) and integrating 
over r' results in 



dt 



2 

/^aVtIa - nxU + (-l)^-^^(/i2 - /Aa/3)VC 

mim2n 



(86) 



Using equation (fTHj) . we identify the right hand side as —V ■ ijixu + jx/rrix) and find 
that 

JA = - DxVpx - (-l)^^^^(/l2 - hap)Vc 

mim2n 

= (-1)^^^ [Dx - xx{h2 - fxap)] Vc. (87) 
mim2n 

where we have used equation (fTTjl to write Vpi at constant pressure and temperature 
(which is the case here) as [p + c{dp/dc)Tp]'^c = p^/(m2n)Vc, and similarly Vp2 = 
— p^/(min)Vc. Since jx = (— l)^p-DVc (by definition), we see from equation (|H7jl that 
the mutual diffusion constant is given by 

D = -^^^[Dx - xxifu - fxap)] = ^/l2. (88) 

1711171211 mim2^ 

where in the last step we used equation for the case A = 1 and equation (jS^ for the 

case A = 2. The Green's function expression for the mutual diffusion in equation (jHHjl 

is therefore identical to the Green-Kubo result in equation ()29|) . 

So far in this section the system was assumed to be in local equilibrium at time 

t = 0. We are really interested in systems that are initially far from equilibrium. If 

left unperturbed, such a system would still reach local equilibrium, after a time that we 
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shall denote by and which depends on the details of the initial condition (with no 
general way to determine it other than through measurement). 

A natural question is therefore whether the Gaussian part of the Green's function 
for such a far from equilibrium system is in agreement with hydrodynamics for times 
t ^ as well. We will now argue that this is indeed the case. 

For this, we start from equation (jHTj) for dnx/dt. That equation was obtained by 
using 1) the Gaussian approximation for the Green's function, and 2) that t is large 
enough for Ka2 to have reached its long time value 2D\. Since the assumption of 
local equilibrium was not needed for the derivation of equation (|HT|l. it holds also in 
cases where the system starts far from equilibrium. To obtain hydrodynamics from 
equation (|HT|l. we calculated K^i assuming that the system was in local equilibrium 
at time 0. Now we would like to do the same in situations where local equilibrium is 
established only later, at time t^. The local equilibrium condition at time ti can be 
used by writing K^i as (cf. equation 



Kai = C'^Ai t Ar' = / dr -— — ((5(r - rxi{tL))) (89) 

j {d{r" - rxi{tL))) 

Assuming now that all microscopic correlations between quantities at time and at later 
times have died out before t^, the condition that t'ai(O) = r' in {vxi(t)6{r" — rxi(ti)))xr' 
should not matter, i.e. 

{vxiimr" - rxi{tL)))xr' ~ {vxiimr" - rAi(tL))) (90) 
Using this in equation (|H^ gives 

Kai = f rfr"(^Ai(t)).,,(*,)=."(5(r" - rxi{tL))) (91) 



where the subscript appended to the average {vxi{t))r^^(tL)=r" indicates that it is taken 
only over all configuration for which rxi{tL) = r" . The quantity {vxi{t))r^^(^tL)=r" can 
be calculated using local equilibrium, since the condition involves time t^, at which time 
we assume local equilibrium to have been established. In fact, except for the time shift 
over ti, it is precisely the same quantity Kai that was calculated for the case where local 
equilibrium was present at time zero, with the result (jHH), so that 

(^Al(t))..,(t,)=r." = U- (-l)^^^(/i2 - fxc.p)^C. (92) 

This expression for (DAi(t))rAi(iL)=r" is independent of r", so that after its substitution 
into equation ()9H). the integration over r" can be performed giving 

kxi = u- - fxap)Vc. (93) 

This is identical to the case where local equilibrium was established at time zero, i.e. 
equation (|S^. Substituting equation in equation (|HT|) would thus lead, as before, to 
the hydrodynamic equation (jHUj) . with the right value for the mutual diffusion constant 
as given by equation (jHHj) . 

In conclusion, the Gaussian approximation to the Green's function for a far from 
equilibrium initial situation gives hydrodynamic behaviour for long times t > ti as long 
as local equilibrium is established at some finite time ti (not necessarily zero). 
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5. Discussion 

In this paper we have shown that for binary diffusion in an isotopic or ideal mixture, 
the Green's function theory is able to handle long time behaviour. 

We remark that the very short time behaviour was discussed in a previous paper 
[20] , using a recent theorem regarding particle displacements and their cumulants under 
the condition that the initial velocity distribution is Gaussian while the initial position 
distribution may be arbitrary. 

The long time behaviour was checked here by expanding in gradients and using a 
local equilibrium assumption for long times. The latter assumption is natural because it 
is generally expected that local equilibrium is a condition for hydrodynamic behaviour. 

Given the form of the Green's functions formula in equation ()36|) . in which the 
density of component A at time t is expressed solely in terms of that same density at 
time zero, it may be surprising that the Green's functions do not produce only self- 
diffusion. In fact the Green's functions reproduce mutual diffusion with the correct 
(Green-Kubo) expression for the mutual diffusion constant (see equation (jSHI)) with 
both the self (Darken-Hartley-Crank) terms and the necessary cross terms fxap and fu- 
Note that the cross terms found in equations ()88j) arose from the local equilibrium 
assumption applied to the drift terms Kxi{r',t), whereas the self-terms come from 
the rather simple approximation to the variance Kx2{'r'',t) ~ 2Dx- We can therefore 
conclude that although the Green's functions seem to treat the system as self-diffusion, 
the sensitivity of the drift terms to local equilibrium restores the true mutual and 
hydrodynamic character of the binary diffusion. 

We note that the Green's function formalism presented here can be straightfor- 
wardly generalized to general classical multi-component mixtures with a non-equilibrium 
initial distribution, since the equations of the Green's functions, equations (IHU jl - lfTSj) . do 
not rely on the fact that there are only two components, nor on the fact that they are 
isotopic or ideal. It was mainly for the discussion of the approach to hydrodynamics in 
section H] that we restricted ourselves to isotopic and ideal binary mixtures in order not 
to make the manipulations overly complicated. 

We will end with a few questions to be addressed in the future. 

Although we expect that the Green's functions will also be able to reproduce 
hydrodynamics at long times in other systems, especially in general mixtures, this is 
at present an open problem. 

The Green's function formalism presented in section El only describes mass 
transport. Analogous to the treatment in references [7,8], it would be desirable to 
extend the formalism presented here to include local currents and energy densities. 

The intermediate time behaviour of the Green's functions for mixtures — expected 
to be more relevant for nanoscale systems — is to be checked in computer simulations. 
Here, the correction terms Bxn in equation ()65|) will become important. If the situation 
is similar to that of the heat pulse studied in references [7-9], taking the Gaussian 
approximation alone would turn out not to suffice for the correct description of the 
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behaviour of the system, but that including just a few correction terms, such as in 
equation (fTSjl . should be enough to get good agreement with the simulations on the 
picosecond scale. Note that these correction terms represent non-Gaussian behaviour 
of the particle displacement; thus it is expected that non-Gaussian behaviour occurs 
on the time scale in between the initial infinitesimal and the long time, hydrodynamic 
behaviour. 

These issues are currently under investigation. 
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A. Local equilibrium 

In this appendix we will formulate the local equilibrium distribution in an expansion 
in gradients following the method of Ernst [27] for a one-component fluid. For a 
one component system Ernst derived an expression for the (grand canonical-like) local 
equilibrium phase space distribution function at time t expanded in gradients of the local 
thermodynamic fields T (temperature), u (fiuid velocity) and fi (chemical potential), 
at a position r and time t at which one wishes to calculate a local quantity. We will 
use a similar derivation here for a two component mixture, in which case there are two 
chemical potentials. As in the main text, we will consider a case with only gradients in 
the chemical potentials fi\, while Vp = 0, VT = and = 0. 

The local equilibrium distribution will be defined in analogy with the usual grand 
canonical equilibrium distribution function at temperature T, chemical potentials /ii 
and fj,2 and an average fiuid velocity n = 0, given by 

^/T.^ _ exp{-[g(r) - ELif^xN,]/ksT} 

where the Hamiltonian is defined in equation Q and the normalization Z = 
En, En, I exp{-[H{r)-fiiNi -fi2N2]/kBT} /{NMlh^''). For local equilibrium, 
we suppose that at a certain instant, taken to be time 0, one can define a uniform fluid 
velocity u and spatial varying chemical potentials of the two components, /^i(t') and 
fj^2{f), while p and T are constant. The local equilibrium distribution is then 

, _ exp{-[g,(r) -Jdr' El=i f^xir')n,{r')]/kBT} 
where the normalization Zx is given by 



Ari=0 Af2=0 
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Furthermore, in equations ()A.2|) and ()A.3|1 the phase functions corresponding to the 
local number densities and the total peculiar energy are (cf. equations ((H) and 

i=l 

2 Nx 2 N„ 

A=l i=l u=l j=l 

We note that the requirement of constant pressure puts, through the Gibbs-Duhem 
relation ()13|) . a restriction on the functions ^Ji\{r): 

V/i2 = -ni/n2V/ii. (A. 6) 

Following Ernst we expand the distribution function as 

VL{T\^lx) = Po(r|^A) + n(r|/iA) + . . . (a.7) 

Here ^\ = fix{r,t) with r a specific position around which local quantities are to be 
evaluated. Furthermore, the leading term Vo(T\fi\) = PL(X\fJ'\)\^^=^^(^rt)y which is a 
distribution of the form of a grand canonical equilibrium distribution in equation ()A.1|) 
at the (fixed) values fix = fi\{r, t) and in terms of peculiar velocities i.e. with H replaced 
hyH,. 

Expanding equation ()A.2|) around fix = fixii^,t), we find for the next-to- leading 
term in equation ()A.7jl 



n(r|/XA 



5VUT\fi> 



d 

[r' - r) ■ Vfixir,t) - t—fix{r,t) 



where we used that fix{r',0) = fix{r,t) + {r' — r) ■ 'Vfix{r,t) — t{d/dt)nxir,t) + . . . . 
The functional derivatives are, using equation ()A.2j) . 



SMTlfix) 
Sfixir') 



VL{r\fix 



nxir 



51nZL 



(A.9) 



kBT 5fix{r')_ 

The term 5 In Zi^jdjixir') in equation ()A.9|) can be shown to be the average of fixir') / ksT 
taken with Pl, so, as it occurs in equation ()A.8|1 . 

5PL(r|/iA_) 



5fix{r') 



knT ^ 



^o(r|/iA).. , M , 



(A.IO) 



where we used that the average ()l ~ ()o is translation invariant so {nx{r'))Q = nx- We 
substitute equation ()A.10|) into equation ()A.8|) and obtain 

^o(ri^A) " ' 



Vi{T\fix 



I dr'Y,[nx{r') - n,][(r' - r) • V fix - t^Ux]- (A.ll) 

A=l 

In equation ()A.11|) . we are only interested in djix/dt to leading order in the gradients. 
This is given by the Euler equations which are just equations ffT^ -(| ^ without 
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dissipation terms i.e. Dp/Dt = Dc/Dt = 0. The quantities fix are just functions of 
p and c (and of T, which is constant), so that also Dpx/Dt = and equation ()A.11|1 
becomes 

Vi{r\px) = / Eb^^Alr) - nx]Vfix ■ (r' - r + ut). (A.12) 

^ A=l 

As in the main text we use instead of pi the quantity p = pi/rrii — p2/fn2. with 
the help of the Gibbs-Duhem relation ()A.6|) . which gives V/i2 = — (ni/n2) V/ii, we find 
VyUi = mi(l — c) V/i and VyU2 = — m2cV/i. Substituting this into equation (jA.12jl and 
using equation (fT^ gives the result used in the main text: 

VoiT\px) 



-PiiTlpx) 



dr'[{\ — c)mini(r') — cm2n2{r')]'V p ■ {r — r + ut) 



dp\ MM f dr' [(1 - c)mini(r') - cm2fi2{r')]{r' -r + ut)-Vc 

n2{r') 



Tp 

Vo{T\px) p I dr' 



ni[r ) 

cm2n 



c)min 



[r' - r + ut] ■ Vc. 



(A.13) 



B. Derivation of equation (I84j) 



The average ()o in equation (jHHj) is over a variant of the grand canonical ensemble where 
each particle has an additional velocity u compared to the normal grand canonical 
equilibrium ensemble. Since each particle has the same velocity shift, this does not 
change the relative motion of the particles. Hence, we can replace ()o in the above 
expression by a real canonical {u = 0) ensemble average (), if we replace vii{t) by 
u + vii(t) in equation (jHBI): 

t) = —^(N^5ir' - rn(0))[it + v^it)] 



niir',0) 



X 



P 



dr" 



ni{r 



n2{r" 



{r" -r + ut) - Vcjy (B.l) 



cm2n (1 — c)min 
The part of the average proportional to u is precisely equal to ni(r', 0), so 



ni(r',0) 
+ p{Ni6{r' -mmvnit) 



dr" 



ni[r") 

071211 



Mr") 

[1 — cjniin. 



(B.2) 
r + ut) ■ Vc 



niir',0) 



To simplify this expression further we will apply the following trick, based on the fact 
that {vn{t)) = 0, several times. The term {Ni6{r' — rii(0))i'ii(t)) can be shown to 
vanish as follows. Its only spatial dependence is through r' but because of translation 
invariance it cannot really depend on r' so we can integrate over r' and divide by V. 
This leaves us with {NiVii(t))/V. Then i>ii(t) can be replaced by i'ii(O) because both 
the number of particles and the grand canonical distribution function are invariant under 
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time evolution and the average of it is then seen to be zero, i.e. {NiVii{0)) /V = 0. Thus 



Mn(r',t) = 1* 



ni 

X 



ni{r") 



(r 



cm2n 



c)min 



iVi5(r'-rn(0))vn(t) 

ut)-Vc, (B.3) 

where we also replaced ?2i(r',0) by rii = ni{r,t), since the difference can be expressed 
in terms of gradients and would contribute only to second order in the gradients in 
equation ()B.3|) . The term proportional to (— r + ut) can be integrated over r". Using 
/ dr"nx{r") = Nx, this yields factors of and A^iiV2 times Vii{t)S{r' - rii(O)). By 
the same reasoning as above this is zero, and we are left with 



Mii(r',t) = u + 



P_ 

ni 



rfr"(iVi5(r'-rn(0))«n(t) 



071271 



n2[r' 



[1 — c)min 



r"-Vc. (B.4) 



The remaining integral over r" can be performed by using the explicit forms of ni{r'' 
and h2{r") in equation ()A.4|) . This gives: 



Mii(r',t) 

= n+^(iVi«n(t)5(r'-rn(0)) 
rii 



N2 



cm2n 



c)min 



) • Vc. 



^ i=2 " i'=l 

Here we have written out the term i = 1 explicitly. This term is in fact zero: once the 
factor of r' has been taken out of the average, we can use the above reasoning that the 
remainder by translation invariance does not depend on r' and can be integrated over 
r', after which we have an integral of a time invariant expression times Vii{t), which is 
zero. Using also that the particles of each component are indistinguishable, we find 



Mii(r',t) =u + 



P_ 

ni 



{6{r' - rn(0))iVi(iVi - l)«n(t)ri2(0)) 



cm2n 

{6{r'-r^,{0))N^N2V^,{t)r2m) 



Vc. 



(B.5) 



(1 — cjmin 

This whole expressions is independent of r'. It may seem to change under translation 
over some fixed vector s as this would change rxi to rxi + s. However, the prefactor in 
front of s that results is itself independent of r', can therefore be integrated over r', and 
is then zero because of the factor Vuit). Hence the expression is translation invariant 
and may be integrated over r' and divided by V. This then gives equation 
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